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Abstract 

I propose that stiffness may be defined and quantified for nonlinear systems using Lyapunov 
exponents, and demonstrate tlie relationship that exists between stiffness and the fractal 
dimension of a strange attractor: that stiff chaos is thin chaos. 



What constitutes a stiff dynamical 
system? Stiffness is closely related to 
numerical methods [1]: the signature 
of stiffness in a problem is that, upon 
integration with a general numerical 
scheme — a method not specially de- 
signed for stiff problems — the rou- 
tine takes extremely small integration 
steps [2], which makes the process 
computationally expensive. One view 
is that stiffness is inextricably linked 
with the numerical integration scheme 
used, so that there would be no such 
thing as an intrinsically stiff dynami- 
cal system, and the best we could hope 
for is an operational definition such as 
that above [3]. Moreover, it has been 
proposed that chaotic problems cannot 
be stiff [4]. I argue below that this is 
not the case, and provide a definition 
and a quantitative measure of stiff- 
ness for nonlinear dynamical systems. 
I demonstrate how stiffness affects the 
geometry of the strange attractor of a 
chaotic system: that stiff chaos is thin- 
ner — has smaller fractional part of 



the fractal dimension 
chaos. 



than nonstiff 



When integrating a stiff problem with 
a variable-step explicit numerical inte- 
gration scheme, the initial step length 
chosen causes the method to be at 
or near numerical instability, which 
leads to a large local truncation error 
estimate. This causes the numerical 
routine to reduce the step length sub- 
stantially, until the principal local trun- 
cation error is brought back within its 
prescribed bound. The routine then in- 
tegrates the problem successfully, but 
uses a far greater number of steps than 
seems reasonable, given the smooth- 
ness of the solution. Because of this, 
computation time and round-off error 
are a problem when using conventional 
numerical integration techniques on 
stiff problems, and special methods 
have been developed for them. 

Traditionally in numerical analysis, a 
linear stiff system of size n is defined 
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by [5] 

Re(Aj) < 0, l<i< n, (1) 

where Aj are the eigenvalues of the Ja- 
cobian of the system, with 



max |Re(Ai)| > min |Re(Ai)|. (2) 

l<i<n l<i<n 



The stiffness ratio R provides a quan- 
titative measure of stiffness: 



R 



max |Re(Aj 

l<i<n 

min |Re(Aj 

l<i<n 



(3) 



By this definition, a stiff problem has 
a stable fixed point with eigenval- 
ues of greatly differing magnitudes; 
large negative eigenvalues correspond 
to fast-decaying transients e^* in the 
solution. 

As an example of a linear stiff problem, 
consider the equation 



y" + 1001?/'+ lOOOy = 0. 



(4) 



We can write this as the vector equation 
y' = Ay where 







1 



-1000 -1001 



(5) 



the e~^°°°* transient term decays, but 
in fact the presence of the large nega- 
tive eigenvalue A2 prevents this. With 
appropriate initial conditions, one can 
even set 5 = and remove the e-^^"*^* 
term from the solution entirely; one 
nevertheless has to use a very small 
step length throughout the calculation, 
as step length is still dictated by the 
size of A2. 

The definition of linear stiffness of 
Eqs. (2) and (3) is not relevant for 
nonlinear systems. The stiffness ratio 
Eq. (3) is often not a good measure 
of stiffness even for linear systems, 
since if the minimum eigenvalue is 
zero, the problem has infinite stiff- 
ness ratio, but may not really be stiff 
at all if the other eigenvalues are of 
moderate size. The inadequacy of the 
stiffness ratio is clearly recognized by 
numerical analysts, who have moved 
away from trying to pin down the def- 
inition of stiffness and have generally 
adopted the pragmatic approach I al- 
luded to earlier: "Stiff equations are 
problems for which explicit methods 
don't work" [3]. 

Let us look at a nonlinear stiff prob- 
lem. In Fig. 1, 1 show the results of in- 
tegrating the van der Pol equation 

y"-fx{l-y^)y' + y^O (7) 



and the eigenvalues are Ai = —1 and 
A2 = —1000. This has solution 

y = Ae-* + Be-^°°°*, (6) 

so when integrating the problem nu- 
merically with a general variable- step 
method we would expect to be able to 
use a large integration step length after 



using a standard variable- step fourth- 
order Runge-Kutta code [6], for the 
cases = 1 and // = 100. There is 
the same bound on the principal lo- 
cal truncation error estimate in both 
cases. We can see from the far greater 
number of steps needed for fi = 100, 
that at large /i the van der Pol equation 
becomes very stiff. The steps are so 
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Fig. 1. Numerical integrations using 
a variable-step fourth-order Runge-Kutta 
method. The principal local truncation er- 
ror is e = in all cases. The van 
der Pol equation with (a) /x = 1 and (b) 
/X = 100; (c) forced van der Pol equation 
with /X = 100, A = 10, and a; = 1. 



small at // = 100 that the individual 
crosses representing each step merge 
to form a continuous broad line on 
the graph. At these large ji values, the 
equation describes a relaxation oscilla- 
tor [7] . These have fast and slow states 
in their cycle which characterizes the 
jerky motion displayed in Fig. 1(b). 
And if we introduce forcing into the 
van der Pol equation 



/-Ml 



v'^)y' + V 



A cos ujt (8) 



we obtain a chaotic system [8-11] 
which, like its unforced counterpart, 
is stiff, as demonstrated in Fig. 1(c)), 
where the steps are so small that they 
have merged together in the image. 

Throughout the foregoing, I have illus- 
trated stiff systems with reference to 
their effect on the numerical method 
used for integrating them. Is it possible 
then to define nonlinear stiffness in- 
dependently of the numerical method 
used? Lambert [5] puts forward the 
definition: A system is said to be stiff 
in a given interval x if in that inter- 
val the neighbouring solution curves 
approach the solution curve at a rate 
which is very large in comparison with 
the rate at which the solution varies 
in that interval. This definition does 
not make reference to the numerical 
method used and is instead concerned 
with the curvature and local Lyapunov 
exponents of the solution. The curva- 
ture of the solution curve 



y" 



[I + |/'2)3/2 



(9) 



at a point quantifies the wiggliness of 
the trajectory at that point, while local 
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Lyapunov exponents are defined as 

-, (10) 



1, ai{t + T} 



7i(T, t) — lim - In , , 



where (Tj(t) are the principal axes of 
an ellipsoidal ball evolving in time in 
phase space. The local Lyapunov expo- 
nent 7i(r, t) is a generalization of the 
Lyapunov exponent A^: the Lyapunov 
exponent is the limit of r going to in- 
finity in the local Lyapunov exponent 



\i = lim ji{T,t). 



(11) 



Lyapunov exponents show the rate 
of convergence or divergence of 
neighbouring trajectories, and an n- 
dimensional system has n Lyapunov 
exponents corresponding to n expand- 
ing or contracting directions in phase 
space. Since Lyapunov exponents are 
defined in the infinite time limit, they 
cannot reflect differing rates of conver- 
gence or divergence in different parts 
of a trajectory. Whereas the Lyapunov 
exponent is the same for almost all 
starting points on a trajectory, the local 
Lyapunov exponent can vary depend- 
ing on the starting point and the length 
of trajectory examined. Fast conver- 
gence to a neighbouring trajectory, 
which implies having large negative 
local Lyapunov exponents, indicates 
stiffness. A system that is equally stiff 
at all points along a trajectory has a 
constant convergence rate, and in this 
case, the local Lyapunov exponent 
will be the same as the Lyapunov ex- 
ponent. Often, though, a system can 
show intervals of stiff and nonstiff be- 
haviour. Relaxation oscillators like the 
van der Pol oscillator are examples of 
this, having a stiff slow manifold, and 



a nonstiff fast manifold. With each 
oscillation we have two intervals of 
stiff behaviour interspersed with two 
intervals of nonstiff behaviour. 

As it stands, Lambert's definition of 
stiffness will not do for us, since it as- 
sumes that the system is not chaotic: 
the principal local Lyapunov exponent 
is large and negative to obtain fast con- 
vergence of neighbouring trajectories. 
Instead, we need to adapt it to allow 
for chaotic behaviour along with stiff- 
ness, by looking at the largest negative 
local Lyapunov exponent, rather than 
the principal local Lyapunov exponent. 
Our definition of nonlinear stiffness is 
then: A system is stiff in a given inter- 
val if in that interval the most nega- 
tive local Lyapunov exponent is large, 
while the curvature of the trajectory is 
small. A quantitative measure of non- 
linear stiffness at any point can then 
be obtained from the ratio of the most 
negative local Lyapunov exponent and 
the curvature of the trajectory: 



R 



■nl 



mm 7i(r,t) 

l<»<n 



K{t) 



(12) 



If desired this could be averaged over 
the trajectory to give a measure of 
mean stiffness. 

We now have a definition of nonlinear 
stiffness; what does it imply for chaotic 
systems? If a Lyapunov exponent of a 
system is large and negative, then the 
local Lyapunov exponent must be large 
and negative at least over some of the 
trajectory, so, given suitable bounds on 
the curvature of the trajectory, a large 
negative Lyapunov exponent is a suf- 
ficient condition for stiffness. On the 



301 



other hand, at least one positive Lya- 
punov exponent is necessary for chaos. 
Thus for stiff chaos we should have a 
large spread of Lyapunov exponents, 
with at least one positive, and one large 
and negative. 

Lyapunov exponents are certainly re- 
lated to the fractal dimension of an at- 
tractor. The Kaplan- Yorke conjecture 
[12] holds that the fractal dimension is 

DKY^TT^j^^k+j, (13) 

k=l 

where j is the expansion dimension of 
the system: the largest integer such that 

E Afe > 0. (14) 

k=l 

The Kaplan- Yorke estimate is often 
good — i.e., very close to the measured 
(capacity) dimension. In general the 
fractal dimension D lies between the 
expansion dimension and the Kaplan- 
Yorke estimate [13]: 

j<D< Dky. (15) 

For example, in a three-dimensional 
chaotic system, Ai > 0, A2 = 0, and 
A3 < 0, so 

Dky ^2 + ^. (16) 

From this, if IA3I gets larger, with in- 
creasing stiffness, then the fractal di- 
mension of the strange attractor moves 
closer to two. 

The size of the fractional part of the 
fractal dimension of a strange attractor 



may be termed the thickness or thin- 
ness of the chaos. The folded structure 
of the foliations of the strange attrac- 
tor — the Abraham and Shaw bagel 
[14] — is tightly wound in the case of 
thin chaos with small fractional part. 
On the other hand, if the fractal fold- 
ing of the surface of the strange attrac- 
tor is macroscopic, the chaos is thick 
chaos having a large fractional part to 
its fractal dimension. 

We should expect then that stiff chaos 
will be thin, and this indeed proves 
to be the case. In Fig. 1 I illustrated 
the stiffness of the forced van der Pol 
equation. This equation is also well 
known for the difficulty of capturing 
pictorially its chaotic nature at large /j, 
[15]. On the other hand, with the Shaw 
variant of the forcing [16], which may 
be written as a three-dimensional au- 
tonomous system 

X— —y + — y^)x 
y — x-\-A cos (jjz 

i = l (17) 

whose divergence — the sum of the 
Lyapunov exponents — is //(I — y"^), 
a strange attractor is easily found at 
A = 0.932, = 1.18, and u = 1.86. 
Now, with the benefit of the above 
analysis, we see how fractal dimen- 
sion, dissipation, and stiffness are all 
linked together: whereas the normal 
forced van der Pol equation studied 
at large pL (high dissipation) is stiff, 
making the chaos thin, the attractor 
in the Shaw variant at small /i (low 
dissipation) displays thick chaos, and, 
following these arguments, will not be 
stiff; this is indeed the case (Fig. 2). 
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Fig. 2. The Shaw version of the forced van 

der Pol equation, parameters A = 0.932, 
/X = 1.18, and lo = 1.86. (a) Poincare 
map showing strange attractor. (b) Nu- 
merical integration using a variable-step 
fourth-order Runge-Kutta method [6], 
with the same principal local truncation 
error e = 10^^ as in Fig. 1. 
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